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Abstract. We study a class of processes that are akin to the Wright-Fisher model, with transition prob- 
Cn abilities weighted in terms of the frequency-dependent fitness of the population types. By considering an 

^~H approximate weak formulation of the discrete problem, we are able to derive a corresponding continuous weak 

C _ ^ formulation for the probability density. Therefore, we obtain a family of partial differential equations (PDE) 

Cn for the evolution of the probability density, and which will be an approximation of the discrete process in the 

K^ joint large population, small time-steps and weak selection limit. If the fitness functions are sufficiently regu- 

_^ lar, we can recast the weak formulation in a more standard formulation, without any boundary conditions, but 

^^ supplemented by a number of conservation laws. The equations in this family can be purely diffusive, purely 

■^^ hyperbolic or of convection-diffusion type, with frequency dependent convection. The particular outcome will 

depend on the assumed scalings. The diffusive equations are of the degenerate type; using a duality approach, 
we also obtain a frequency dependent version of the Kimura equation without any further assumptions. We 
also show that the convective approximation is related to the replicator dynamics and provide some estimate 
of how accurate is the convective approximation, with respect to the convective-diffusion approximation. In 
particular, we show that the mode, but not the expected value, of the probability distribution is modelled by the 
replicator dynamics. Some numerical simulations that illustrate the results are also presented. 
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Evolution is naturally a multiscale phenomenon (Keller, 1999; Metz, 2011). The choice of right scale 
to describe a particular problem has as much art as science. For some populations (e.g, with non overlap- 
's^ ping generations) a discrete time provides adequate description; for different examples, this is excessively 
■^ simplifying. Large populations can be described as infinite (in order to use differential equations, for ex- 
•/^ ample), but this imposes limitations in the time validity of the model (Chalub and Souza, 2009b). On the 
' , other hand, some finite population effects, like, for example, the bottleneck effect, will be missing in any 
t^^ description relying in infinite populations (Hartle and Clark, 2007). 

^^ In this vein, diffusion approximations, frequently used for large populations and long time scales, enjoy 

, a long tradition in population genetics. This tradition dates back as early as the work by Feller (1951) and 

. . references there in. In particular, diffusion approximations were implicitly used in the pioneering works of 

Wright (1938, 1937) and Fisher (1922, 1930). These efforts have been further developed in a number of 

directions as, for instance, in the studies on multispecies models in Sato (1976, 1983); see also the review 

%^ in Sato (1978). Subsequently, Ethier and Kurtz (1986) systematically studied the approximation of finite 

^ Markov chain models by diffusions. In particular, they showed the validity of a diffusion approximation 

to a multidimensional Wright-Fisher model, in the regime of weak selection, and linear fitness. This led 

to a notable progress in diffusion theory, as reported for instance in (Ethier and Kurtz, 1986; Stroock and 

Varandhan, 1997). This considerable progress, in turn, led to a large use of diffusion theory in population 

genetics, which can be verified in contemporary introductions to the subject (see Ewens, 2004; Etheridge, 

2011). 

There is also a more heuristic approach, called the Kramers-Moyal expansion, where the kernel of the 
master equation of the stochastic process is fully expanded in a series. The diffusion approximation can 
be viewed as a Kramers-Moyal expansion truncated at second order. Although it is commonly claimed 
that the full expansion is needed in order to obtain a continuous approximation of discrete processes, it is 
known that under various conditions discrete Markov chains can be approximated by diffusions; cf. Ethier 
and Kurtz (1986) and Stroock and Varandhan (1997) for instance. In this work, we shall show that under a 
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number of conditions similar results hold for the discrete processes considered. See (Van Kampen, 2001) 
for a discussion about this and other techniques for continuous approximations of discrete processes. 

As observed above, results along similar lines had been obtained earlier by a number of authors (Feller, 
1951; Ethier and Kurtz, 1986; Ewens, 2004). These works approach the problem mostly within a prob- 
abilistic framework, while here we take a pure analytical setting, and this brings two immediate conse- 
quences: firstly, we are able to derive a weak formulation for the forward Kolmogorov equation, assuming 
only continuity of the fitness functions. The second, and possibly most important consequence, is that we 
are able to deal with a variety of scalings for the evolution problem. This yields a full family of evolution 
problems: genetic-drift dominated evolution, which is described by a diffusion equation; selection domi- 
nated evolution, which is governed by a hyperbolic equation; and an evolutionary dynamics, where the two 
forces are balanced, which is governed by a convection-diffusion equation that we term replicator-diffusion. 

If we assume some more regularity of the fitness functions, we can then recast the weak formulation in a 
strong formulation. In this case, we cannot impose any boundary conditions, but we must supplement it by 
a number of conservation laws, namely that the probability of fixation of each type for a given probability 
density of population, in any time, must be the same as in the initial time. The conservation laws are used 
to circumvent the impossibility of imposing boundary conditions when the boundaries are absorbing. 

Furthermore, by a duality argument we obtain the backward equation formulation. For the particular 
case of linear fitness and balanced scalings, we then recover the classical result by Ethier and Kurtz (1986). 
Additionally, by an appropriate combination of the weak and strong formulations, we are able to give a 
complete description of the forward solution. 

A complimentary approach to the study of evolution, based on evolutionary game theory, has also been 
developed (cf. Smith, 1982) with conclusions that are not always compatible with results from diffusion 
theory. As an example, diffusion models without mutation lead to the fixation of a homogeneous popula- 
tion, while frequency dependent models associated to the replicator dynamics' may lead to stable mixed 
populations. For an introduction to evolutionary game theory and replicator dynamics, we refer the reader 
to Hofbauer and Sigmund (1998) and Weibull (1995). 

Consistent interaction among these two modelling schools have been attempted by a number of authors, 
with different degrees of success (see Traulsen et al., 2005; Lessard and Ladret, 2007; Lessard, 2005; McK- 
ane and Waxman, 2007; Waxman, 201 1; Champagnat et al., 2006, 2008; Fournier and Meleard, 2004). We 
will show, as in many of these works, that both descriptions — the one based on the diffusion approxima- 
tion and the one based on the replicator dynamics — are correct as models for the evolutionary dynamics 
of a given trait, but in different scalings. As a byproduct, we will provide a generalisation of the Kimura 
equation valid for an arbitrary number of types and general fitnesses; we will suggest that the replicator 
equation is a model with limited time validity, given a certain maximum admissible error, and that the 
solution of the replicator equation indicates the most probable state (mode) to find a population, not the 
expected value of the trait^. 

The work presented here is a development of earlier work in (Chalub and Souza, 2009b, a): the for- 
mer studying the derivation and convergence of the Moran model with two types to the 1 -d version of the 
replicator-diffusion equation discussed here, and the latter with a comprehensive analytical study of the 
1-d replicator-diffusion equation. The derivation of the continuous model in Chalub and Souza (2009b) 
hinged on the idea that a formal expansion of master equation, but with control of the local error, and 
results on well-posedness of the continuous classical problem can be brought together via numerical anal- 
ysis approximation results. This combination then yielded uniform convergence, in any proper closed sub 
interval of [0,1], of the rescaled probabilities of the discrete model to the continuous probability density. 
This convergence result, combined with the analytical results in Chalub and Souza (2009a) on a weak for- 
mulation that satisfies the conservation laws provided a continuous measure solution. The discrete process 
then converges weakly towards such a solution, on a neighbourhood of each endpoint, but uniformly as 
described above. To study the Wright-Fisher continuous limits, however, we took a different route. We 
split the neutral and selection parts of the discrete process. This allows us to derive an approximate discrete 
weak formulation of the discrete process, with global error control. Further, by embedding the discrete 



In this work, we will use the expressions "replicator dynamics", "replicator equation" and "replicator system" indistinctly. 
We call "mode" the most probable state in intS""', and not in S""' (see definition 2). 
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Detailed model 


Meaning of parameter 7 


Simplified model 


Kinetic models 


mean free path 


hydrodynamical models 


Othmer-Dumbar-Alt model 


mean free path 


Keller-Segel model 


Quantum Mechanics 


rescaled Planck constant 


Classical Mechanics 


Relativistic mechanics 


(rescaled light velocity)" ' 


Non-relativistic Mechanics 


Moran process 


inverse of population size 


replicator-diffusion equation 


Moran process 


inverse of population size 


replicator equation 



Table 1 . Detailed and simplified models. The last two lines state that both the repli- 
cator equation and the replicator diffusion equation approximates the Moran process. 
References to these works are (Bardos et al., 1991, 1993; Cercignani, 2002; Hillen and 
Othmer, 2000; Othmer and Hillen, 2002; Chalub et al., 2004; Stevens, 2000; Hepp, 1974; 
Cirincione and Chernoff, 1981; Bjorken and Drell, 1964; Chalub and Souza, 2009a,b). 



probabilities in an appropriate measure space, we could use compactness arguments to obtain the continu- 
ous limit. Thus, in this setting both the weak formulation and the weak convergence of the discrete model 
to the continuous one follows with considerable less effort, but we do not get the improved convergence on 
the interior. 

1.1. Scalings, limits and approximations. In order to be able to study somewhat more general models, 
we follow the approach used by the authors in Chalub and Souza (2009b). In particular, we are interested 
not only in diffusion approximations, but in approximations that can be consistent with the dynamics of the 
corresponding discrete process. 
We begin with a definition: 

Definition 1. We shall say that a simplified model ^q is an approximation of the family of detailed models 
^y, 7 > 0, in a sense %, where % is an appropriate metric as, for instance, any norm in a suitable space of 
functions (e.g., L^, L^, L°°, etc) if the following holds true: 

(1) Consider a certain family of initial conditions hL such that limy_^o/Zy = h^, in the sense %; 

(2) Evolve through the model ^y the initial condition hi, and through the model ^o the initial condi- 
tion /zq until the time t < °o obtaining hy{t) and /zo(f ) respectively; 

If for allt <°° we have that limy_^o hy{t) — ho{t), in the sense %, then we say that the model ^y converge in 
the limit 7— > 0, in the sense %, to the model ^q. If, furthermore, this convergence is uniformly in t G [0,oo), 
then we say hat the model ^y converge in the limit y ^f^to the model ^0 uniformly in time. 

Some examples of the relation between detailed and simplified models are listed in Table 1 . 

In general, some extra assumptions are frequently required to allow the passage to the limit. If, for ex- 
ample, there are more than one small parameter in the detailed model, it is natural to assume a relationship 
among them, called scaling, as, in general, the limit model will depend on how these parameters approaches 
zero. Other assumptions may also be necessary, as it will be discussed in the next paragraph. The process 
of taking the limit of a family of models, considering a given scaling, will be called "the thermodynamical 
limit"; by extension, we shall also call the limit model the thermodynamical limit. In this work, depending 
on the precise choice of the scaling, the limit equation can be of drift-type (a partial differential equation 
fully equivalent to the replicator equation or system), of purely diffusion type, or, in a delicate balance, of 
drift-diffusion type. 

In what follows, an important and natural assumption that must be introduced in order that we have an 
approximation in the sense of definition 1 is the so-called weak selection principle, to be precisely stated 
in equation (11). Generally speaking, we assume that the fitness of a given individual (i.e., the probability 
of finding descendants of this individual in the next generation) decreases to 1 when the time separation 
between two successive generations At approaches zero. This is a natural assumption when we consider 
that two successive generations collapses into a single one. However, in most of the literature, the weak 
selection principle is assumed in the limit of A^ ^ 00, where A^ is the population size. Although they are 
equivalent (as we shall assume a certain scaling relation between A^ and At), we consider our approach 
more natural. 
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In this work, we will consider as the detailed model, the Wright-Fisher process, to be rigorously intro- 
duced and analysed for finite populations in section 3: an evolutionary process for an asexual population of 
A^ individuals, constant in size, divided in n different types, that evolves according to a specific rule, with 
fixed time separation between generations of A/ > (the detailed model in the discussion above, where 
7 is the inverse of the population size — or, as we shall see, equivalently, the inter generation time). In 
short, given a certain scaling and the weak selection principle, we find a certain partial-differential equation 
of drift-diffusion type with degenerated coefficients (what we call the replicator-diffusion equation) as the 
thermodynamical limit of the Wright-Fisher process. However, if we are interested only in the first time 
scale of the Wright-Fisher process, we shall assume different scalings and obtain as simplified limit the 
replicator equation, a first order ordinary differential system. 

We prove that the limit f ^ oo of both models are well defined and limy_^o lim,^oo hyit) — lim^^oo limy_^o hy{t). 
We therefore, conjecture that this convergence is uniform in time for t e IR+. 

1 .2. Outline. Section 2 introduces the basic notation and provides an extended abstract of our main results. 
In section 3, we review some classical results about the discrete process (the finite population Wright- 
Fisher process); we also show the existence of a number of associated conservation laws, and obtain an 
asymptotic representation for the second order moments, in the large population regime. In section 4, we 
obtain a family of continuous limits of the Wright-Fisher process depending on the scalings that are derived 
within a weak formulation, with solutions in appropriate measure spaces. In particular, we derive the 
replicator-diffusion equation, and show that it satisfies continuous counterparts of the conservation laws for 
the discrete process. We then continue the study of the replicator-diffusion equation in section 5, where we 
derive the main properties of its solutions, including a description of the solution structure as a regular part 
and a sum of singular measures over the sumbsimplices, and the large time convergence to a sum of Dirac 
measures over the vertices of the simplex. We also show that the probability distribution associated with all 
types in the population concentrates along the evolutionary stable states. Additionally, in subsection 5.2, we 
obtain the backward equation as the proper dual of the replicator-diffusion equation, providing a consistent 
generalisation of the Kimura equation for the n types and arbitrary fitness functions. In section 6, we study 
the replicator equation and show that, in the regime of strong selection the solutions to the replicator- 
diffusion will be well approximated by the solutions to the replicator equation within a finite time interval. 
Numerical examples are given in Section 7, where we also point out that, for intermediate times and large 
but finite populations, the replicator equation will approximate the mode of the discrete evolution, but not 
the expected value of a given trait. Conclusions are presented in section 8. 

2. Preliminaries and main results 
We begin by introducing the space of states for the evolution: 
Definition 2. Let M+ ~ [0,°°). We define the n — 1 dimensional simplex 



5"-^=<^xeR'! 



!=1 



We also define the set of vertices of the simplex A5"^ := {x G 5"^ |3i,x,- — 1}, its interior int5"^' :— 
{x € 5"^ jV/jX; > 0} and its boundary dS"^ — 5"^ \int5"^'. The state of the population is a vector 
X G S"^ . The elements ofAS"^^ are denoted e,-, i ~ !,...,« and called "homogeneous states". A vector 
X € 5"^ \A5"^ is a "mixed state". 

In what follows, we let p{x,t) to be the probability density of finding the population at state x e S"^^ at 
time f > 0. 

Definition 3. We call the fitness of a given type a continuous function Xj/^'' ; 5"^ — > M, and the average 
fitness in a given population is given by V^(x) :— Y4=iXiW i"^)- Note that we consider the fitness function 
to not depend explicitly on time. 



n 

Strong selection in this context is not directly related or opposed to weak selection as introduced before. 
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In this work, we derive a family of detailed models described by a parabolic equation of drift-diffusion 
type, with degenerated coefficients (DiBenedetto, 1993; Carrol and Schowalter, 1976), defined in the sim- 
plex S"^\ called the replicator-diffusion equation, namely: 

r d,p - i^„_i,,/, := f rijii dij {D,jp) - r,zl dt {^,p) , 

(1) } Dij-.^Xidij-XiXj , 

\ a:=x,-(v/«(x)-V/(x)) , 

with / , j = 1 , . . . , n — 1 , K" > 0, and where 5,/ = 1 if / = y and otherwise is the Kronecker delta. The above 
equation has a solution in the classical sense (i.e., everywhere differentiable). Furthermore, in the classical 
sense, it is a well posed problem, without any boundary conditions. However, this classical solution is not 
the correct limit of the discrete process. In order to find the correct limit, equation (1) is to be supplemented 
with n conservation laws. From now on, whenever we refer to the replicator-diffusion equation (1), we are 
implicitly assuming these conservation laws. 
Our main conclusions are: 

(1) An analysis of the equation (1) leads to a unique solution of measure type. This will require 
definitions of appropriate functional spaces. 

(2) This unique solution approximates, in the thermodynamical limit, the evolution of a discrete popu- 
lation by the Wright-Fisher process pointwise for any time. In addition, the large time asymptotics 
is consistent with the discrete model. 

(3) A reduced model, obtained by setting K" = in (1) (with only one conservation law), is shown to be 
equivalent to the replicator dynamics. This will suggest that the replicator dynamics approximates 
the discrete process for any f, however with an error increasing in f along a fixed discretisation 

(4) Furthermore, the solution of the replicator equation models the time evolution of the mode of the 
probability distribution associated to the discrete process (and not the expected value of the same 
distribution); 

(5) A frequency dependent generalisation of the Kimura equation for an arbitrary number of types is 
obtained by looking at the dual problem for (1). 

Before going into the technical details, we explain the last paragraph a little further 
Equation (1) has two natural time scales, one for the natural selection (the mathematical drift and, as 
we shall see, fully compatible with the replicator equation), the second for the genetic drift (the mathemat- 
ical diffusion). That is why we call equation (1) together with the conservation laws to be introduced in 
subsection 4.5, the "replicator-diffusion equation". More precisely, the solution of the replicator-diffusion 
equation when K ~Q (which is of hyperbolic type) is the leading order term of the solution /?«: of the 
replicator-diffusion equation for small K (i.e., large fitness and/or short times). The replicator-diffusion 
equation with zero diffusion (k" = 0) happens to be the replicator equation (or system) (Hofbauer and Sig- 

mund, 1998). In an appropriate sense, to be made precise in section 6.3 (theorem 11), we have p^ — > po, 
pointwise, but not uniformly in time. 

This theorem cannot be made uniform in time, for general fitness functions and initial conditions, as the 
Wright-Fisher process always converge in f ^ oo to a linear combination of homogeneous states, while it 
is possible that the solution of the replicator equation converges to a stable mixed state. 

The former is the mathematical formulation of a known principle in evolutionary biology that states 
that "given enough time every mutant gene will be fixed or extinct." (Kimura, 1962). This means that 
the final state of the replicator-diffusion equation with any K" > will be a linear combination of Dirac 
deltas at the vertices of the simplex 5"^'. Actually, for any positive time, the solution of equation (1) 
with the conservation laws described above is a sum of a classical function in the simplex plus a sum of 
singular measures over all the subsimplexes on dS"^^ and, inductively, on their boundaries subsimplexes. 
In particular, we shall have also Dirac measures supported on the vertices of the simplex. These measures 
appears immediately, i.e., for any f > 0. This represents the fact that in a single step there is a non zero — 
however, small — probability that the population reaches fixation through Wright-Fisher evolution. The 
full evolution and the final states of the replicator-diffusion equation will be studied in section 5. 

From the practical point of view, we are, however, often interested in transient states ("in the long run, 
we are all dead", said John Maynard Keynes), specially because the transient states become more and more 
important for the discrete evolution as the population size increases. Heuristically, when the population is 
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Figure 1 . The boxes in the figure represents the solutions of three different models: the 
Wright-Fisher process (finite population A^), the replicator-diffusion equation (positive 
diffusion k) and the replicator equation. The vertical axis indicates the arrow of time 
(top-down), and the horizontal axis indicates, first the large population limit, secondly 
the no-diffusion limit. Consider that there is a maximum error acceptable error e (in the 
L°° norm) between the Wright-Fisher model (suitably Radonmised — see subsection 4.3) 
and the continuous approximation. Therefore, there is a population size A'o such that for 
N > No the difference between the replicator-diffusion equation and the discrete model 
is less that e. For the replicator approximation, for any A^, it may exist maximal time 
fmax(^) < °° such that for t > fmax that the error is too large. 

large the stochastic fluctuations decreases in importance, and therefore, its evolution is deterministic. The 
associated limit will be given by equation (1), with K" = 0, i.e., the hyperbolic limit of equation (1). This 
equation does not develop finite-time singularities. This is one more peculiarity of equation (1): diffusion 
seems to be a deregularizing effect; the solution of the parabolic replicator-diffusion equation K > is 
less regular than the solution of the hyperbolic null-diffusion limit K" = 0, contrary to most examples in 
the literature (John, 1991; FoUand, 1995). See, on the other hand, (Murray, 2003), for diffusion driven 
instability. 

The relationships between the three models is summarised in Figure 1. 

Finally, we observe that the natural formulation for the continuous limit is the weak one. For such a 
formulation, we only require the fitness functions to be continuous. If, in addition, these functions are 
also Lipschitz, we can then recast the problem in a strong sense, provided that we supplement it with the 
conservation laws. Finally, requiring the fitness functions to be smooth allows for a number of results 
about the solutions to be easily derived. In particular, one can prove a structure theorem that show that 
the problem is equivalent to a hirearchy of classical degenerate problems, provided that some members are 
interpreted as densities for singular measures. 



3. The discrete model 

In this section, we study the discrete model, i.e., the Wright-Fisher model for constant population, arbi- 
trary number of types and arbitrary fitnesses functions. We start, in subsection 3.1 with basic definitions; 
in subsection 3.2 we briefly review some important results in the literature. We also prove that the discrete 
process has as many conservation laws as types. Additionally, we also show that the final state is a linear 
superposition of these independent stationary states, with coefficients that depend on the initial condition 
and that can be calculated from a set of n linearly independent conservation laws. All these results will be 
useful in the correct determination of the continuous process, to be done in sections 4 and 5. The discrete 
Wright-Fisher process was studied, with different level of details in, for example, (Ewens, 2004; Nowak, 



THE FREQUENCY-DEPENDENT WRIGHT-FISHER MODEL: DIFFUSIVE AND NON-DIFFUSIVE APPROXIMATIONS. 7 

2006; Imhof and Nowak, 2006), but, to the best of our knowledge the conservation laws associated to the 
process were overlooked. 

The fact that the final state in the Wright-Fisher process, among others, for a finite population is always 
homogeneous was also a matter of dispute with respect to the validity of the modelling (Vickery, 1988; 
Smith, 1988). As we will shortly see in this work, this dispute is basically a consequence of the existence 
of two different time scales hidden in the model: the non-diffusive (drift) and the diffusive one. See 
also Ethier and Kurtz (1986) and Etheridge (2011) and references therein for a discussion in the role of 
time scales. 

3. 1. Preliminaries. We consider a fixed size population of A^ individuals at time t consisting of a fraction 

X, e {0, ^ , ^ , . . . , 1 } of individuals of type i — 1 , 2, • • • , n. The population evolves in discrete generations 
with time-step separation of Af. We introduce the following notation: 

Definition 4. The state of a population is defined by a vector in the N -discrete n — l-dimensional simplex 

S"f^-^ ■.= lx={xi,--- ,x„)\\x\:='£xi = l,Xie |o,-,-,---,l| I . 



We also define the set of vertices of the n — l-dimensional simplex 

an- 



A5«-i := {x e Sl-'\3i,xi = 1} = {e,|/ =!,...,«} 



The elements of AS"^^ are called the homogeneous states. To each type we attribute a function, called 
fitness, ^^ : S'^ — > (0, °°). It is convenient to assume that ^^^ is a discretization of a smooth function on 
the simplex S"^ ; more assumptions on ^^ will be introduced in section 4. 

A population at time f + Af is obtained from the population at time t choosing A^ individuals with proba- 
bility proportional to the fitness. More precisely, we define the average fitness ^^(x) ~ ^f^jX,^^^ (x) and 
then the transition probability from a population at state y to a population at state x is given by 

{Nxi ) ! (A^X2) ! • • • (A^x„) ! fj^ y ^'^ (y) J 

The evolutionary process given by a Markov chain with transition probabilities given by equation (2) is 
called the (frequency dependent) Wright-Fisher process. 

Let ^(f) = (P(x,f))^g_5„-i, with 

^eT:={P:S"^-'xR+^R+\ ^ P(x,-) = 1}, 

where P{x,t) is the probability of finding the population at a given state x G S"^ at time t. Then, the 
evolution is given by 

(3) P(x,f+Af) = (^^(f))(x):- Y. @NMy^x)P{y,t). 



3.2. Stationary states, final states and conservation laws. We call an homogeneous population a popu- 

nn- 



lation of a single type, i.e., P{x,t) ^ P^{x) for v G AS^j ', where 



A(y) 

From the inner product definition: 



1 , y = x 

0, y^x 



{v,w):^ E ^(x)w(x) 

XGS" 



s^' 



it follows immediately that {Px,Py) = dx,y = 1 if x = y and otherwise. 

Now, we state classical results for the Wright-Fisher process that will be useful in the sequel. The 
interested reader should consult Karlin and Taylor (1975). 
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Lemma 1. A function f defined in S"^ is a fixed state of the operator ,9' if and only if f is a linear 
combination of homogeneous states. In particular, S' has exactly n linearly independent eigenfunctions 
associated to the eigenvalue X = 1. For all non-negative initial condition P^, the final result is a linear 
combination of homogeneous states, 

n 

P~:=limP(.,f) = E^i[^']4- 

Definition 5. We define a linear conservation law as one given by a linear functional L over the functions 
of S"^ such that \-{^{t-\- At)) ~ L(^(f)). A set of linear conservation laws is linearly independent, if 
the only linear combination providing a trivial conservation law L(^(f)) = is L = 0. 

Proposition 1. Define F''' :~ Lxes"-' Fp Px, ' = 1, • • • ,n, a functional over S^f . Therefore F(''(x) is the 

fixation probability associated to the initial condition x. Finally, the set {f^',..., f^"'} is a basis for the 
set of linear conservation laws associated to the operator ^. 

Proof. Given any initial condition P e T, we define Fp as the fixation probability of the type / in a popu- 
lation initially in the state P. From the fact that 



5^~P=£FJ')pe, 



i=l 



we find 



In particular 



Finally, 



^W ^ (^oop) (g.) ^ (^~p,4) ^ (p, (^^P,^) . 



L 4''Px= E {p.,{,^ype,)p.. 



{^T^= E ^/'^x = F('). 



Therefore, F*^'^ is an eigenvector of ,3^^. In particular, 

F«(e;) = {{3rYP,,,P,^) = (Pe,,=^%) = {Pe.,Pe,) = 5, 



V ■ 



It is immediate to prove that they are linearly independent; let ai , . . . , a„ such that Y!i=i o;,F('' = 0, i.e., for 
every x e ^J^T^ L"=i «iF(')(x) = 0. Using x = e,-, we conclude that a,- == 0, and then {F^^', . . . , F^"'} is a 
basis for the eigenspace of ^^ associated to A = 1. 

Now, consider a linear conservation law L. From standard representation theorems, there is a vector 
we 5j^"' such that 

(4) {^{t),w) = L(^(f)) = L(^(f +Af)) = (^^(0,w) = (^(0, .3^M . 

Therefore, w is an eigenvector of ,3^^ associated to A = 1 and then it is a linear combination of F*^'', 
i — !,...,«. 

D 

Remark 1. The conservation of probability (the most natural conservation law), follows directly from the 
equation 

EF»(x) = E^/' = l'Vxe5«-^ 



THE FREQUENCY-DEPENDENT WRIGHT-FISHER MODEL: DIFFUSIVE AND NON-DIFFUSIVE APPROXIMATIONS. 9 

3.3. Properties of the transition kernel. The probability conservation is a consequence of the defini- 
tion (2) and reads 

(5) Y. ®NMy^x) = l, VyG5«-i. 

It also follows from the definition that 

which can be readily interpreted by the absence of mutations in the model. 
We also observe that we can write: 

(7) ©^^Ar(y^x)=Afy,x,A^-z)s(y,x,A^,Af), 



where 



i. A'! A NX, 



(9) srv.x.A?.Af):=ni ^^^^ 1 



S(y,x,A^,Af):=n -^^ 



Remark 2. The factorisation (7) splits the transition kernel in two factors: the fitness dependent E and the 
fitness independent A (more precisely, in the neutral case ^^ independent of i, we have that S = 1, and 
therefore — A.}. As it will be seen in the next sections, the continuous approximation will naturally present 
two time-scales, the first one associated to the drift (in the mathematical sense) - or natural selection 
- whose precise expression will depend on the asymptotic expression of SA (see lemma 3); the second 
time-scale is the diffusive one - the genetic drift - and will depend solely on the first three momenta of 
A. Therefore, on one hand, the splitting in "natural selection" and "genetic drift" is already present 
in the transition kernel, i.e., in finite populations; on the other hand, any process such that both fitness- 
dependent and fitness-independent parts have the same asymptotics will have the same thermodynamical 
limit, possibly in different scalings. This might explain why the replicator-diffusion equation (1) seems to 
be a robust approximation of many different evolutionary processes (Ewens, 2004; Ethier and Kurtz, 1986; 
Etheridge, 2011; Der et al, 2011). 



It will be also convenient to write 
and to introduce 



S"i^^^{yeR"-'\x±yeS%-'}. 



1 " 

zTi^yt, z^—^. and ^x.r = {t e M"| V T; = and|Ti| <x;/z}. 

Lemma 2. For large N (and then small z) the neutral transition probability A has the following first 
moments for X G int5"^^." 



^ A(x,x + zT,z) = l , 
^ T,-A(x,x + zT,z)=0 



N.x 



N.x^ 



^ T;T/A(X,X + ZT,Z) = (x;5y -X;X,-) + 0{z^) . 






where 5ij ~ 1 ifi^j and otherwise is the Kronecker delta. 

Proof. The first one is a simple consequence of normalisation. The second one indicates that the average 
displacement in the neutral case is identically zero. 



10 
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For the third one, we use the StirUng formula xl — \/2nxx^e ^ (l + 0{x ')) to write 



A^! 



{27ty 



N- 



{Nx,)\{Nx2V.---{NXn)l N"-^ ixiX2---XnM''xl'''---xf'' 

We then write 



1+0 



A^ 



Since the right hand side is smooth, standard estimation of Riemann sums yields 

^^ E T.T,A(X-ZT,X,Z)=/ T,T,n(^y"dT[l+0(z2)] 
[271) 2 ,^^^.„-l J.yx,z •■-■ V Xi 



TV, A ■ 

•^x,z ■ \ 1^=1 



(=1 



Xk-ZTk\ , __ ,_l Xn+zZlJiT^k 



■y^,. 



T,T,exp<^ - ^ 



m + 2 



Xk 



n-l ^m+2 



-x„log 



i-i rn+l ' ;ji+l I 2- "^^ 



m+2 



A:=l'^<r 



a=l 



|jjdT[l+0(z2)] 

dT[l + 0(z2)] 



T,T/exp< -, 



«— 1 fl.2 1 tt— 1 ( n—\ ^ 



m=l / 



onde P„i{x) is a linear combination of monomials in T;, / = 1 , . . . , n of degrees 3 to m. From the symmetry 
of the integrand and the domain (in the odd case), we find 

f , if « odd , 

T,T,A(x-zT,x,z)T„,T„,...T„„dT=<; ^^ ^J,^n(.-\^^\ if„even. ' 



■y^..i 



G,„+2(l+0(e-T^ 



where Gm+2 G K+, and (J„_i is the smallest eigenvalues of the quadratic form ^ to be shortly defined. 
Therefore, 

TiTyA(x-ZT,X,z)dT 

(l+0(z2))dT + 0(e-^ 



TiTyA(x-zT,X,z)dT= / T;TjA(x-zT,X,z)dT 

(27r 



I -II 



^Xl ...x„ jk"-! 






-It? 1 "r' '""' 



= ^^^^ / T,T,exp|-i^(T,T)ldT + 0(z2) , 

i2 is the quadratic form associated to the matrix matrix F = {Fij), i,j = 1, • • • ,n — 1, defined by Fa 



x^ and Fij —x^, for ij^j; its eigenvalues are denoted by ai > a2 > • • • > (7„_ i > and CJi • • • (7„_ 



(xi---x„) 



-1 



/ T;T,A(x,X + ZT,z)dT=-'-^|^ / T;T,-e ^^("'^'dT 



Finally, we use that 
(10) 



^J^?^372 / T,T,( tr.n-/ )e-^^f^^^W + 0(z2). 
2iUUxi)'^^ -M"-^ '\ti i/k J 



XiXjQ Z-^(^''')dT = [2%) V ^xiX2...X„(x;5,7 - XiXj) 



[ T,-T,T4e-Z'^('''^'dT = , V/, j,/t == 1, . . . ,n 
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and finish the proof. For the detailed calculation of equation (10), see appendix A. D 

D 



4. Continuations of the discrete model 

The aim of this section is to obtain a differential equation that approximates the discrete evolution, when 
the population is large (A^ -^ °o) and there is no time-separation between successive generations (Af -^ 0). 
The relevant variables, x e 5"^' and f > will be forced to be continuous. 

The first four subsections will be devoted to the development of three models, based on partial differen- 
tial equations obtained from the Wright-Fisher process, when N -^ °° and Af — )• (see equations (19), (20) 
and (!'), respectively). There is no "right choice" of the simplified model. As we could expect, simpler 
models will have a restricted application. For example, the model given by equation (19) is equivalent to 
a system of a ordinary differential equations; actually, it is exactly equivalent to the well-know replicator 
dynamics (see Hofbauer and Sigmund, 1998). On the other hand, the diffusive approximation, given by 
equation (20), is a parabolic partial differential equation that is much simpler to solve than the full model; 
in fact, explicit solutions are know using Gegenbauer polynomials (Ewens, 2004). Our focus will be on the 
replicator-diffusion approximation, equation (1'), which we expect to be valid uniformly in time. 

Results known for the Wright-Fisher process, and stated in section 3 will guide the derivation, i.e., 
the choice of the right thermodynamical limit. We start in subsection 4.1 by the asymptotic expansion 
of the transition kernel in the negligible parameters (suitable combinations of A^ and Af); we plug this 
expansion into the master equation (3) in subsection 4.2. In subsection 4.3, we construct the continuous 
version of the discrete probability densities; in particular, we interpolate discrete probabilities in order to 
represent them by continuous probability measures; these measures will be central when we finally pass 
to the limits in subsection 4.4, obtaining the various continuous approximations of the discrete model. 
Finally, in subsection 4.5, we show that, for every conservation law of the discrete process, there exists a 
corresponding conservation law in the continuous model. As a by product, the final state of the continuous 
model shall be a linear superposition of homogeneous states (see lemma 1 and compare it with theorem 7). 

4.1. Preliminaries. From a biological point of view, the most important assumption in this derivation is 
the so called weak selection principle 

(11) ^«(y) = l + (AO>»(y), 

where \j/^'' : S"^ ' ^ M is a continuous function, and v > is a parameter yet to be specified^. In this case, 
we also have 

*Ar(y) = l + (Af)Xy). 

Next, we obtain asymptotic information about the sums of the transition probabilities of the whole 
process. This will turn out to be a key result in the derivation of a continuous limit. 



Lemma 3. For large N and small Af, with \/N{AtY small, with A and E given by equations (8), (9) 
and (11), it is true that 

^ S(x,x + zT,A^,Af)A(x,x + zT,z) 






(=1 



l+zN{AtyY,Ti(xi/'^''>{x)-xff{x))+0{zN{Atf\{At)^\N\At 



A(x,x + zT,z) 



See also the discussion on the weak selection principle and the choice of At (and not N) in its expansion in Chalub and Souza 



(20()9b). 
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Proof. Initially, using lemma 2, we find that 



(12) 



L L ^i^J (v^'''(x) - r(x)) (/'^(x) - V/(x)) A(X,X + ZT,Z) 






: ^ (x,5„ -x,x,) (v/«(x) - v/(x)) (/'^(x) - v/(x)) + 0{z^) 
:fx,(v/»(x)2-v/(x)2)+0(z2). 



Using an asymptotic expansion in the equation (9), we find 

S(x,x + zT,A^,Af) 



V/«(x)-v/(x) 



: expJA^f (x,-+zT,){(Af)^ [v/«(x) - v/(x)] - (Af)^^ 
:exp|A^f (x,-+zT,){(Af)^ [v/«(x)-v^(x)^ ^^'^'' 



(Af 



,2v 



V/»(x)v/(x)-V/(x)2]+O((A0^^)} 



V.(-)(,)^(,)_^(,)2+ (V^^"W-r(x))^ 



V/«(x)2-v/(x)2)+0((Af)3^)} 



-0((Af)3v)} 



: expi zA^ (Af )^ f T,- (va(') (x) - v/(x) j - ^ 



A^(Af 



2v Af 



^Xi h/W(x)2-v/(x)2 +0(zA?(Af)2^A^(Af)3^) 



1 +zN{Aty t ^' (V^'H^) r(x)) - ^^^ Ex, (v/(')(x)2 - v/(x) 



1=1 



1=1 



2^r2/'A,^2v n 



Z^N^iAt) 



ij=l 



E T,T4r»(x)-V/(x) U(-'''(x)-v/(x) +0(zA^(Af)2^A^(Af)3^A^2(^^ 



Finally, we multiply by A(x,x + zT,z), add in zT G 5" + and using equation (12), lemma 2 and z^A^ = 1, 
we finish the proof. D 

D 

4.2. An asymptotic weak-discrete formulation. For any A^, we have the master equation: 
p{x,t+At,N)= E 0A',A((x-y^x)/5(x-y,f,A^). 

■' N,\- 

Let g : S"^^ X M — > M be an appropriate test function. On multiplying by g{x,t) and summing over S'Jf , 
we find that: 

E p{x,t+At,N)g{x,t)= E E 0A',AKx-y^x)p(x-y,f,Af)g(x,f) 

E E 0A',A((x^x + y)/5(x,f,A^)^(x + y,f) 

E /'(x,f,A^) E A(x,x + y,z)S(x,x + y,z^Af)^(x + y,f) 
e^r' yes'"-' 

E Pi^^t^N) E A(x,x + zT,z)S(x,x + zT,z^Af)g(x + zT,f) 



on— 1 



xGs;, 



2TG5", 



We can now use our asymptotic information on the process as follows: 
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Proposition 2. Let g be the restriction of a C ' (i2), where £1 is any open set such that 5"^ C ii. Then, 
under the same hypothesis of Lemma 3, we have 

(13) Y. {p{^,t + At,N)-pix,t,N))g{x,t) 



xes", 



L /'(x,f,A^) 



XGS" 



1 «— 1 n—l 

^ Xi{5ij^Xj)dljg{x,t) - (AtY ^ xjd,jg(x,t){xj/^^'\x) - v/(x)) 



2N 



ij=l 



.7=1 



+ 0{zN{Atf'',NiAt)^'',N\Aty\z\z^iAty) 
Proof. For sufficiently large A^ (and then small z), we may write: 
^ p{x,t+At,N)gix,t)= Y, p{x,t,N) 



XGS" 



xes". 



ZTGS 



l+zNiAtyY^i[w'''H^)-W{^)]+0{zN{Atf\NiAt)^\N\Atf'') 



'n,x+ 



A(x,x + zT,z) 



«-I 



^2 n-I 



g(x,f)+Z^T,-5^.g(x,f) + - ^ TkTldl,^g{x,t)+Z^R{x,T,t,z) 



j=l ^ kj=l 

where there exists a constant C such that: 



\R{x,T,t,z)\<C{l + \\Tf). 



Hence 



Y pix,t+At,N)gix,t)^ Y /'(x,f,A?)< L A(x,x + zT,z)^(x,f) 



xes", 



XGV, 



.zres", 



-z Y A(x,x + zT,z) 



ZTGS" 



i=l 



A^ (Ar)^ Y r<'' (x) - V/(x) T, + ^ T,5,^^(x, f 



n-1 



.7=1 



-z^ Y A(x,x + zT,z) 



Z-fGS" 



n-1 

E 



H n—l 



E ^^irt^(x,0 + E E ^(AO'^.,^(x,r)(V/«(x) - V/(X))T,T; 

'=ii=i 



-z^ ^ A(x,x + zT,z)/?(x,T,f,z) ^+0(zA^(Af)2^,A?(Af)3^A?2(^j 



z^esl 



We now analyse the coefficients of the z power expansion. 
For the leading order, we have 



Y A(x,x + zT,z)g(x,f)=§(x,f). 



zxeS 



N,x+ 



By virtue of Lemma 2, the first order coefficient is zero, and third order coefficient is (9(1). We now 
compute the second order coefficient, which we divide in two terms: the first one is 



n-1 



Ti;T/ 



Y A(x,x + zT,z) Y -Y^v,g{^.t) 



ztesi 



kj=\ 



n—l n—l 

Y{xk{l-Xk)dlg{x,t))- Y {^kxidl^,g{x,t)) 

k=l k,l=l,k^l 



+ 0{z'). 
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The second one is 

^ A(x,x + zT,z)£^A^(Af)^5,^^(x,f)(v/«(x)-v/(x))T,T,- 






= N{Aty 
= N{Aty 



n n—\ 

IE 

f "f ^,^.g(x,f)(v/«(x) - V/(x))(x,5„ -x,x,) + 0{z^) 

'X:'x,5.,^(x,0(/'^(X) - V/(X)) - E-^/^-;^(X'OE-^Kr''Hx) - V/(X)) +0(z2 

;=1 7=1 i=l 



fi-1 



= A^ (Af )" ^ x,-5,^^(x, f ) ( v/(-''' (x) - v/(x)) + O (zA^ (Af ) 

7=1 

Combining the calculations above, we obtain the result. 



D 
D 



Equation (13) can be seen as discrete weak formulation for p{x,t,N) in space only, and thus any lim- 
iting argument would require some regularity assumption on p{x,t,N) in f. In order to circumvent such 
assumptions, we need a full discrete weak formulation: 

Proposition 3. Let T ~ MAt, where M is some fixed positive integer, and let g be an admissible test 
function, with support in S"^ x [0, T]. Let 



J^{kAt}, k = Q,...,M-l. 



Then we have that 



(14) Y. Y. pix,t,N)igix,t+At)-gix,t))- Y Pi^,t,N)gix,0) 



'eTxGsr' 



XGS 



■n-l 

'N 



E E /'(^'f'^) 



'eTxGsr' 



1 

2N 



n-\ 



Y Xi{Sij-Xj)d^jg{x,t) - (Aty Y Xjd^jg{x,t){\j/^J\x) - \j/(x)) 



\2v-l Ar^A*N3v-l »r2/ 



7=1 

\-l _2/»,NV-h 



+ 0{zN{Aty''-\N{Aty''-\N\Atr-\z\At)-\zHAty-'). 

Proof. Sum (13) over T, and estimate the error term by its total sum, taking into account that there are 
(9((Af)^') terms in this sum. This shows the right hand side of (14). To obtain the left hand side, we 
perform a summation by parts and use that ^(x, T) = 0. See appendix B for details. D 

D 

4.3. Continuous representation. The aim is now to obtain a continuous version of (14), but without 
taking any limits yet. We first need some preliminary definitions: 

Definition 6 (Piecewise time interpolation). Let T be a set of sampling times as above, and let To be a 
set of times such that for each f G T, there exists a unique <^ G Tq such that t, G (f,f + Af). Let g be an 
admissible test function with support in S"^^ x [0,r]. Observe that under the assumptions on the sets T 
and Tq, for each t € [0, T] there exists a unique f G T such that t £ [t,t+At), and a unique t, G (f,f + Af). 
With this in mind, we define: 



and 



g(x,t)^g(x,t), fG[f,f+Af), fGT, 



}{x,t)^g{x,£,), te[t,t+At), <^G(f,f+Af), t eT and ^ eTo 



Remarks. For fixedx, we have on one hand that g{x,t) is just freezing the value of g on [f,f+Af) to be 
the value of g{x,t). On the other hand, g{x,t) is freezing the value of g on the same interval to be the value 
of gix^t,), with t, G (f,f+A?). The natural choice for E, willarise, in the present context, from applications 
of the mean value theorem to g over the interval [/,/ + A/]. 



THE FREQUENCY-DEPENDENT WRIGHT-FISHER MODEL: DIFFUSIVE AND NON-DIFFUSIVE APPROXIMATIONS. 15 

Definition 7 (Radonmisation (sic) of discrete densities). Let p{x,t,N) be a probability density defined no 
S'J^ X T. Let 8x denote the atomic measure at x. We define 

PN{x,t)= ^ p{x,t,N)d^, te[t,t+At). 

With this definitions we have the following result 

Proposition 4. Let g be an admissible test function, let N^ ' = K (Af )^, and let p^t (x, t) = P^-i(^)-fi (x, t). 
Then there exists a set To as in Definition 6, such that 

PAt{^,t)dtg{x,t)dxdt~ / /7Ar(x,0)|(x,0)dxdf 

^^/?A((x,f) I Y, ^ii^ij - Xj)dfj8i^^f) 1 dxdf 



K-(Af) 



JS' 



(15) 



(Af 



.v-I 



Js"-^ 



/'A/(x,f) 



n-1 



.i=l 



^x,(/')(x)-v/(x))a,|(x,f 



dx 



+ 0((Af)2^-t-i,(Af)3^-''-i,(Af)T-i,(Af)^+^-i) 



Proof. For the right hand side, we observe that g{x,t) — g{x,t) for x G 5]^^' and IcAt < t < {k+ l)Af, 
k^O,l,..., and that this also holds for all partial derivatives of g not involving t . On using the definition 
of pn, we readily obtain the equivalence between the sums over Sl^ and the integrals in x. For the time 
integrals, we point out that both p„{x,t), g{x,t) and similarly for the derivatives of g are piecewise constant 
in f. Hence the summation over time can be exactly converted into a time integral with a factor of (Af)^'. 
For the error term, a crude estimate of summing de (9((A/)^') terms yield the result. As for the left hand 
side, apply the mean value theorem to g(x, •) to get the result and the set Tq. D 

D 

Remark 4. The reader is cautioned that, although (15) has a remarkable resemblance with a weak for- 
mulation, it is not quite so, since the prospective test functions g and g are not test functions in the usual 
sense. 

4.4. Passage to the limit. We now deal with the limit Af ^ in (15). 

Theorem 1. Under the same assumptions of Theorem (4), we have that, for any choice of parameters jU 
and V, there exists p G L°°{[0,T],BM^{S"^^)) such that pA/(x,f) — > p(x,/) weakly as At — > 0. Moreover, 
the following limits also hold: 

p(x,t)d,g{x,t)dxdt 

n-1 



-JS' 



Js"-^ 

PAti-X-j) 



PAj{x,t)d,g{x,t)dxdt 
Js"-^ 

/'Ar(x,0 f Y. Xi(Ai-Xj)^u8i^>0 1 dxdf 

n-1 , , 

E -^> ( r'-'^ (x) - r(x) ) dig{x, t) dxdt 



Js"-' 



JS' 



p(x,t) Y. MSij-Xj)d^ig{x,t) dxdt 



-JS' 



./'(x,?) 



Yxj(xif^^\x)-xlf{x))dig{x,t) 



dxdt 



Proof. From the tightness of Radon measures, cf Billingsley (1999), we have that, for afixedf,/?Af(x,f) -^ 
p{x,t), lis At -^ 0. 

The convergence of the x-integrals follows from the weak convergence of /?A( -^ p, and from the fact 
that for a continuous function h, we have 



lim ||/z — /zl 



lim ||/i — /li 

Ar^O 



0. 



The convergence of the time integrals follows from the dominated convergence theorem. 



D 
D 
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If either jU < 1 or v < 1, we can multiply (15) by (Af)"™'"'^"'''^"''. It is then easily verified that the 
error term vanishes in the limit, as well as the term with a time derivative. Thus, in this case, we obtain 
stationary limits governed by the steady version of the equations derived below. Now let us assume that 
jU,V> 1. Additionally, we have that \/}V(Af)^ <C 1, as in lemma 3, implies V >jU/2. Therefore, if jU e [1,2) 
or 2v > /i > 2 it is easily verified that the error term will be small. If both jU, V > 1, we have stationary 
solutions given by the initial condition. 

The other cases are as follows: 

Theorem 2. There exists peL°°{ [0, T] ; BM+ (5"" ' ) ) such that 
: If pL E (1,2), V = 1, the convective or drift approximation: 



(16) 



(17) 



(18) 



p{x,t)dtg{x,t)dxdt- / /?(x,fo)g(x,fo)dx 

J5"-l JS"-^ 



n(x,f) 



n-1 



Jo Js"-^ 
If pL = \, V > I, the diffusive approximation 

p(x,t)d,g{x,t)dxdt 
fi-i 



Y^XjUj\x)-xlf{x))djg{x,t) 

7=1 



dx. 



Js"-' 

K 



s"- 



p{x,tQ)g{x,tQ)dx 



2 Jo Js"-^ 



p(x,f) Y. ^ii^ij-^j)^f}gi^i^) dxdf. 



\i,j=^ 



If ^ = 1, V = 1, the case where there is a maximal balance of selection and genetic drift; we find 
the Replicator-diffusion equation 



Js"-^ 

K 



2 Jo J5"-' 



p{x,t)dtg{x,t)dxdt- / /5(x,fo)g(x,fo)dx 

JS"-1 

p{x,t) Y. xi{5ii-xj)dljg{x,t) dxdt 



+ 11 ,P(^,t) 
Jo Js"-^ 



n-1 



Y^jiw^'H^)-Wi^))djg{x,t) 
J=i 



dx. 



Proof. The result follows from theorem 1, and from straightforward bookeeping of the At orders of the 

terms in (15). D 

D 



Equations (16), (17) and (18) are written in the weak form. In population dynamics, and in others 
contexts as well, they are used casted into the strong formulation (or standard PDE formulation) as follows 
(see, however, remark 5): 

• If /i e (1,2) and V = 1, the convective of drift approximation: 



n-1 



(19) 



d,p = ^Y. ^i ""i ( V^*'' (x) - ^^^) ) P 



(20) 



i=i 

This equation is equivalent to the replicator dynamics, showing that the Wright-Fisher process will 
be equivalent to the the replicator dynamics, in the limit of large population and small time-steps, 
if the population increases faster than the time-step decreases. 
• If /x = 1 and V > 1, the diffusive approximation 

K "^^ 
dtP^i^ E ^'/ {(xi5ij-XiXj)p) , 



which is relevant when the fitness converges to 1 as Af — > faster than A^ — > oo. 
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• When there is a perfect balance between population size and time step, i.e., jU = v = 1, we find the 
replicator-diffusion approximation, given by equation (1), which we repeat here for convenience: 

i,j=l 1=1 

We shall focus on the last equation and on its weak formulation (18). 

Remark 5. We shall see below that the weak and the PDE formulations are not equivalent, and that the 
more appropriate formulation is actually the weak one. 

4.5. Conservation laws from the discrete process. Let us write 6 for the set of all functions g : S"^^ x 
[0, +°°) such that there exist an open set Q. D S"^ ' and a function G : i2 x [0, +°°) — > M such that g is the 
restriction of G to 5""' and G G C^' (i2). 

Notice that in the right hand side of (18), p is multiplied by: 

(21) J XI A74^+£a5,g = 0. 

i..7=I !=1 

Equation (21) is readily seen to be a steady backward equation. We now show that the weak solutions have 
also conservation laws. 

Theorem 3. Let p be a solution to (18) (we shall take (17) as a special case). Let (p Cz & be in the kernel 
of (21). Then 



p{x,t)(p{x)dx= p(x,0)(p(x)dx, 

for almost every t G [0, °°). 

Proof Let 77(f) e Q.([0,°c)), with 77(0) = 1. Then 

g(x,t) = ri{t)(p{x) 
is an admissible test function. On substituting in (18), we find that 

p(x,t)(p(x)ri'(t)dxdt+ / p(x,0)(p{x)dxdt ^0. 
Js"-^ Jo Js"-^ 

Since 77 is an arbitrary function with compact support in [0,oo), the result follows. D 

D 

A similar argument shows also the following 
Theorem 4. Let p be a solution to (16). Then 

p{x,t)dx— / /?(x,0)dx, 

for almost every t G [0, °°). 

Therefore, the conservation laws given by equation (4) now become 

(22) -/ p(t,x)(p(x)dx^Q, 

at Js"-^ 

where (p satisfies (21). In principle the condition set out by (22) seems to imply an infinite (likely to 
be uncountable) number of conservation laws. The following result shows that it is actually much more 
conspicuous: 

Theorems. Lef e,- denote the vertices ofS"^. Then there exist unique Pi, i— 1,. . . ,n^l, with Pi{ej) — 5ij 
that are solutions to (21). In addition, let po = 1. Then, any solution to (21) in & can be written as a linear 
combination of pi, i — Q,.. . ,n—\. In particular its kernel, for solutions in 6, has dimension n. 
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Proof. Given a vertex e,, let e^ be an adjacent vertex. Now we solve (21) in the segment ejc; with boundary 
values 5ij. In the segments not adjacent to e, define the solution to be zero. This defines the solution in 
all one-dimensional simplices. For each two-dimensional subsimplex, we now solve the Dirichlet problem 
with the data from the previous step. Now, assume that we have the solution uniquely defined in all 
subsimplices of dimension m. Repeating the construction above yields the solution in all subsimplices 
of dimension m+1. Proceeding inductively, this yields a solution in S"^^ that is unique and admissible. 
Uniqueness follows from the maximum principle applied at each subsimplex level. Let (p G &, and let 

n 

1=1 

Then <i> (^&. By the preceding argument, <I> and (p must agree at all edges of S"^^. Proceeding inductively 
once again yields that ^ = <I> in 5"^ ^ D 

D 

Thus, the solutions to (18) must satisfy: 

(23) - [ pi{x)p{x,t)dx = 0, /=!, •••,«. 

at Js"-^ 

From a probabilistic viewpoint, the p,, /= 1, ...,«— 1, are naturally the fixation probability of type /. 

We now give a pure analytical argument for this fact. In Section 5, we shall prove Theorem 7 which shows 

that the final state is given by 

n 

P"[P^] = lim/7(-,f) = ^;r,[p']5e, , 
i=\ 

where 5e, is a Dirac measure supported on the vertex e, e -SJ^T • 

5e,.(x)0(x)dx = (^(e;) . 

Clearly, ni[p^] is the fixation probability of type / in a population initially described by a probability distri- 
bution p^. 
Therefore, 

n[^] = I p,(x)/5"(x)dx = / p;(x)/?'(x)dx = / Pi(x)^o(x)dx = p,(xo). 

Remark 6. In the neutral case, i.e., Xj/'^'' (x) = \j/^-'' (x) for all i,j = \,...,n and x e 5"^ ', we define the 
neutral fixation probability nf [5^] — x,-, which follows from the fact that in the neutral case, Pi{x) = x,-. 

5. The Replicator-Diffusion approximation 

We now discuss the nature of solutions p to (!') together with the conservation laws (23). The main 
result of this section is theorem 7. This must be understood as the continuous counterpart of the lemma 1. 
We do not refer to the discrete model to prove this result. Our approach is based solely in the properties 
of the partial differential equation (!'), the restriction of the domain to the domain of interest, and the 
associated conservation laws (23). 

An outline of the proof of theorem 7 is as follows: First, we show that a solution to (18) can be written 
as regular part plus a singular measure over the boundary. Moreover, the regular part vanishes for large 
time. Repeating these arguments over the lower dimensional subsimplices, and using the projection result 
in Proposition 5, we arrive at a representation of /? as a sum of its classical solution and a sum of singular 
measures that are uniformly supported on the descending chain of subsimplices of 5"^ ' down to the zeroth 
dimension. Since the solutions over the sumbsimplices also have a regular part that vanishes, we can 
show that all measures that are not atomicly supported at the vertices should vanish for large time. Thus, 
conservation of probability implies that the steady state of (18) is a sum of deltas. 

Finally, we provide two applications. In subsection 5.2, we study the dual equation. This will be the 
continuous limit of the evolution by the dual equation (backward equation) of the discrete process and 
therefore its solution /(k,f ) gives the fixation probability at time f of a given type (to be prescribed by the 
boundary conditions in the dual process) for a population initially at state k. This gives a generalization 
for an arbitrary number of types and for arbitrary fitnesses of the celebrated Kimura equation with reversed 
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time (Kimura, 1962). In the sequel, subsection 5.3, we will show that if one type dominates all other types 
then, for any initial condition, the fixation probability of this type will be larger than the neutral fixation 
probability. This shows, in particular, that for large populations, the most probable type to fixated will 
be the one playing the Nash-equilibrium strategy of the game (assuming the identity between fitness and 
pay-offs, which is standard in this framework). This is not true in general for small populations (Nowak, 
2006). 

5.1. Solution of the replicator-diffusion equation. We now study in more detail the features of the so- 
lution to (18) and show two important results: first that in the interior of the simplex, the solution must 
satisfy (1') in the classical sense; second, no classical solution to (1') can satisfy the conservation laws. 
Throughout this section, we shall have the further assumption that the fitnesses are smooth. 

Lemma 4. Let p be a solution to (18). Let K C 5"^ be a proper compact subset. Then, in K, p satisfies 
(1') in the classical sense. In particular, p is C°"{K). 

Proof. Let g e C^(K), we have then the standard weak formulation of (1') in K. On the other hand, (1') is 
uniformly parabolic in any proper subset. Hence the weak and strong formulations coincide — c.f. (Evans, 

n 

D 



2010; Taylor, 1996). 

We now obtain some more information about this solution on S"^^. 
Lemma 5. Let p be a classical solution to (1 ') over the interior of S"' 

lim/5(x,f) =0, xeK. 

Proof. We write the drift part as i2 = i V0 + b, define jUs(x) = x\X2 



Then 



with /is 
find 



(24) 



if and only if x e ^5" . In the new variable u — jUse 



x„ (such that /is(x) > in 5""^) 
p and after some manipulations, we 



dtu - 



CO 



n-l 

i=i 



n-l 



CO 



Y^DijdjU-Xi b' 



Iv. 


CO 


''-DVu 


CO 




V2 



Bm 




"-{x^bJ 



with CO = e'^'/jUs and B,- = x,- (b' - L^I 

We shall now study the eigenvalue problem associated to (24) by considering the dual problem with 
regularised coefficients: 

1 



(25) 



V- 



CO'' 



-D^'^'Vcp^ 



-ia)(^)B-V(p(^)==A('^)a)(^)(p 



(<:) 



(p^'> = in dS" 



-I 



D uniformly in x, and Og 



v^r 



where d('^'(x) is a positive defined matrix in S" ^ with D^^' — 
with jUg > in S"^^ and jx^ — > jXs uniformly in x. 

(£] 

For any e > 0, the dominant eigenvalue Xq is real and from the maximum principle it follows that 
Xq 7^ 0. For s = 0, Xq is negative and therefore from its continuity in s, we conclude that Xq <0 for 
any s. Therefore, for any other eigenvalue ReA*^*^' < (see (Evans, 2010) for further details). 

Moreover, let e^^ —> be a decreasing sequence of positive numbers, and (p^^'<^ be the normalised eigen- 
functions for the corresponding leading eigenvalues. Since the coefficients are assumed smooth, the eigen- 
functions are also smooth. Hence, by Rellich theorem, there is a subsequence ek such that (p ''J converges 
in L^{S"^^). By considering the weak formulation for (25), we immediately see that, for this subsequence, 

we must also have X *J — > A. Hence, the leading eingenvalue of (25), with e = 0, is also negative, and all 
other eigenvalues will have smaller real parts than the dominant eigenvalue. 
This also shows that there exists a > 0, such that 



1 



d, 



u-'codV 



s"- 



V- 



CO 



1 



-DVm-Bm 



udV < -a 



u'-codV. 
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Therefore 



p^e '''/isdjc^ u^codx'^0 , 



a 

D 



Lemma 6. Equation (24) has a unique solution u G C(5" ) n C°° (int5" ) 

Proof. Let D^^' and jUg as in Lemma 5. For e > 0, (24) is uniformly parabolic, and hence it has a unique 
solution with the required regularity. 
We write (24) in weak form as 



u^^'{t,x)dt(p{t,x)(0^^'AxAt+ 
Js"-^ 

I ©(''M -DVm('')-Bm(''M •V0(f,x)dxdf+ / M('')(O,x)0(O,x)djcdf =0. 
Js"-^ V2 / A"-' 

Notice that m'*^' e Wq' (5"^'), hence by Rellich Theorem (Evans, 2010), we can select e^ ^ such 
that m'*^*' -^ u^*' in L^(5"^^). For a fixed initial condition, the last integral is independent of e. For the 
first integral, since 5"^' is bounded, we have that L^ (5"^') C L'(5"^'); hence m'*^*' converges inL'(5"^^). 
Thus the first integral converges, by the monotone convergence theorem. The remaining integral can be 
seen to converge by compactness of Sobolev inclusions. 

Finally, the maximum principle can be used to show uniqueness, and regularity follows from analogous 
arguments. D 

D 

This last theorem has an important consequence 

Corollary 1. No solution to (1 ') in the classical sense can satisfy the required conservation laws. 

Proof. Since (!') is uniformly parabolic for any proper set of 5"^ ' , it is possible to show that/? e C(5"^'). 
Given e > 0, choose K C 5"^' such that the m(5"^^\/r) < e, m being the Lebesgue measure in 5"^^. Also, 
let /?" e BM+(5"^^) and y/be a conservation law. Then 

/?''(x)v/(x)dx = a >0, 



and 

/ p(x,f)v/(x)dx= / /5(x,f)v/(x)dx+ / /5(x,f)v/(x)dx 

Js"-^ Js"-^\K Jk 

<Ce+ p{x,t)\j/{x)dx<Ce, 
Jk 

for sufficient large f . Hence there is no conservation. D 

D 

Remark 7. As an extension of Lemma 4.1 in Chalub and Souza (2009a), we observe that if p is a Radon 
measure in 5"^\ we can write p = q + r, with sing supp(g') G dS"^^ and sing supp(r) G int5"^'. 

Proposition 5 (Face Projections). Let « > 1 and let p be a solution to (18) and let S"^^ be a face ofS"~^. 
Assume that sing supp(/?) n5"^ 7^ 0. Then, over S" , p satisfies (18) in one less dimension with forcing 
given by the regular part of p evaluated at x, = 0, for a certain value of i. 

Proof. Assume, without loss of generality, that i = \. In view of remark 7, we can write p = q + r, 
where sing supp(^) C S"^^ with the singular support of r lying in the complement with respect to the full 
simplex. Moreover, we can also assume, without loss of generality, that S"^^ is given by the intersection 
of the hyperplane xi = with S"^K Let us write x = (x2, . . . ,x„). Let h be an appropriate test function in 
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S" ^, satisfying /z(x,0) — and let rj{xi) G Cc{[0, 1]), with 77(0) = 1. Then g = rjh is an approprite test 
function for 5"^' and a direct computation with (18) then yields 

p{xi,x,t)dtg(xi,x,t)dxdt 

p{xi,x,t) Y, Xi{5ij-Xj)d^jg{xi,x,t) dxdt 



2 Jo J5"-' 



p{xux,t) 



'n-l 



70 JS"-^ 

Over xi = 0, on using the definition of g, we find that 

K 



^ Xj ( v/*'^ (xi , x) - v/(xi , x) ) <9^'^(xi , X, f ) 
,7=1 



dx. 



n-I 



q{x,t)d,h{x,t)dxdt ^ '-^ I / ?(x,f) V x,(5i; -x,)5,^/z(x,f) dxd? 

75"-2 Z Jo JS"-2 I -^^ J I 



fl(x,f) 



n-I 



^x,,-(v/<^)(x)-v/(x))a,/z(x,f) 

.7=2 



dxdf. 



For r, we have 



Js"-^ 

K 



r{xi,x,t)dtg{xi,x,t)dxdt 



2. Jo 



n-l 



,^_/{xux,t) I Y, Xi{5ij-Xj)d^jg{xux,t) j dxd? 



+ / / r{xux,t) 
JO Js"-^ 



n-I 



dxdf. 



^x,-(v/<')(xi,x)-v/(xi,x))^/^(xi,x,f) 
..7=1 

By Lemma 4 r is smooth. Therefore, the above equation can be integrated by parts to yield, an integral 
on 5"^ ^ that will cancel out identically, since r is a classical solution to ( 1 ' ), and a number of integrals over 
the various faces of S"^ ^ . In particular, at xi = 0, we find that 

2 Jo Js"-^ 

By collecting together the two calculations on xi = 0, we obtain the result. D 

D 







r(0,x,f)/!(x,f)dxdf. 



In what follows, we shall need some preliminaries. Recall — see Stanley (1996) — that to the simplex 
S"^^ is associated a corresponding /-vector, such that the entry / + 1 (/i+i) is the number of /-dimensional 
subsimplices of S"^^. We shall assume that, for each dimension /, there is a definite order of the sub- 
simplices S''-', with / = 0, . . . ,n — 1 and y = 1, . . . ,/i+i. For / > 0, we shall indicate by S''^ the interior 
S''^ . Moreover, we define the adjacent operator by ad( j,k) which denotes the kth adjacent subsimplex of 
dimension / + 1 to S'' . Notice that there are n — i such simplexes. 

Theorem 6 (Solution Structure). Let p be a solution to (18), and let 5'' be the Radon measure with unit 
mass uniformly supported on S'^. Then the solution p can be written as 

(26) pit,x)=p„i+ Y P,l^"^ ^ = Ull^{i}x {!,...,/,}, 

where pij satisfies (18) on Sij, with forcing given by 

k 
where rij is the regular part ofpij. 

Proof. Start with p defined on S"^^ and let p„i be the classical solution guranteed by Lemma 6. By con- 
sidering (18) with test functions with compact support on int5"^\ we see that p — p„i vanishes. Therefore 
sing supp{p — p„i) C (95"^'. By Remark 7, we can write p = Pni+I, with sing supp((j') C dS"^^. Now 
dS"^^ is the union of /„_i copies of S"^^. By Proposition 5, q must satisfies (18) in one less dimension 
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in each of the subsimpUces with the appropriate forcing. Proceeding inductively, we can now choose a 
subsimplex S"^^ of S"^^. Now, we repeat the argument above for each simplex S"^^ which has S"^^ as a 
subsimplex. Iterate until arrive at the simplices of zero dimension to get the result. D 

This theorem leads to the following result: 

Theorem 7 (Final State). Let 

/'°°(x) = lim/?(x,f) , 

where p is the solution of equation (18) (the weak version of equation (1)) subject to conservation laws (23). 
Then p°° is a linear combination of point masses at the vertices ofS"^,i.e, 

(27) p-^f^7ti[p']5e,. 

i=l 

Proof First, we observe that Lemma 5 still holds if applied to nonhomegenous version of (24), provided 
that the forcing decays for large times. The results now follows from a straightforward application of 
Proposition 5 together with Lemma 5 applied in a descending chain of simplices down to dimension 1 . 
Conservation of probability then yields that /?°° must be a sum of atomic measures at the vertices of S"^ ' . 
On using the other conservation laws, we obtain the coefficients, and hence the result. D 

D 

5.2. Duality and the Kimura equation. The formal adjoint of equation (1) (changing the flow of time 
from forward to backward) provides a generalization of the celebrated Kimura equation (Kimura, 1962), 
both including more types and allowing frequency dependent fitness: 



K 



n—\ «— 1 

32. 



(28) d,f = ^l,J := - ^ A74/+ E a<5,/ . 

In diffusion theory this equation is associated with a martingale problem for the diffusive continuous pro- 
cess. In genetics, the meaning of equation (28) is seldom made clear and depends on the boundary con- 
ditions imposed. One possible and common interpretation is as follows: given an homogeneous state 
e, e A5"^', let /i(k,f) be the probability that given a population initially in a well-defined state k e 5"^' 
(i.e., /?'(x) :— /?(x,0) = ^(x)) we find the population fixed at the homogeneous state e,- at time t (or 
before), i.e., /i(k,f) is the probability of having p{x,t) = 5e,(x) from the initial condition 5k(x), given by 
(p(-,f ), 5e,). In this case, we need to find consistent boundary conditions. See Maruyama (1977); Etheridge 
(2011). 

This follows from the interpretation of the discrete adjoint evolution: 

QiKt+At)^ E 0^^A,(k^k')e(k',f), 

which reads as follows: "the fixation probability after a time interval f + Af of a given type, for a population 
initially at state k is equal to the sum over all possible states k' of the fixation probability after a time 
interval f of a population initially at state k' times the transition probability from k to k'". 

Let us study the fixation of type 1 , represented by the state e 1 . Let us now call V, the face of the simplex 
with Xi = (type / is absent). Then, fly =0. For i =^ 1' /i|y. i^ the solution of dtf — -S?,_2 i^f^ where the 
type / was omitted from the equation. As the faces of the simplex are invariant under the adjoint evolution 
(one more fact to be attributed to lack of mutations in the model), this represent the same problem in one 
dimension less. We continue this procedure until we find the evolution in the edge from vertex 1 to vertex 
! 7^ 1, Lii. In this case, we have that /L : [0, 1] — > M, the restriction of f to this edge, with k the fraction 

\L[j 

of type 1 individuals, is the solution of 

(29) dj='^k{l-k)dlf + k{l-k) (v/<^'L,,« - r*''L,,«) */ 

with boundary conditions given by /(O) =0 and /( 1 ) = 1 and y/"' I is the restriction of !//('' to the edge 
Lii. The forward and backward versions of Equation (29) are fully studied in the references (Chalub and 
Souza, 2009a,b). For Vf^^'l, — V^''' I , constant this is the Kimura equation. 

\Lij \Lij 
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5.3. Strategy dominance. Let us assume that y/'^' (x) > y/''^ (x) for all x e 5""' . This happens, for exam- 
ple, if we identify fitness functions with pay-offs in game theory, types with strategists, and if strategist 1 
plays the Nash-equilibrium strategy. 
Therefore, we prove 

Theorem 8. If, for all states x G 5"^^ and all types i = !,...,«, '^f^ \x) > ^/'^''(x), then the fixation 
probability of the first type is not less than the neutral fixation probability for any initial condition p^; i.e, 

ni[p']>7:f[p']. 

Proof First note that it is enough to prove that 7ii[5x] > 7if[8x] = x\ for all x G 5"^^ The difference 
Pi(x) — xi satisfy 

\ E A7^,7(pi(x)-xi) + £a,5,(pi(x)-xi) = -ni = -xi(v/(^nx)-V/(x)) <0, 
/,/=i (=1 

with vertex conditions p\ (e,) — xi (e,) = for / = !,...,«. Now, we proceed by induction in n. For the case 
n = 2, the proof is in (Chalub and Souza, 2009b, section 4.3); we reproduce it here only for completeness. 
We write explicitly the equation for pi : 

|x(l -x)dlpx +x (^(''(x) - v/(x)) d^px = 

with pi(0) = and pi(l) = 1. We simpHfy the equation using the fact that v/(')(x) - v/(x) = (1 - 
x) ( y/''' (x) — Xj/ (x) ) and the solution is given by 



piW = 



/o^exp[-f/o*(v/(')(l)-V/(2)(l))dl 



dx 



./o exp [- i /o' ( V^(i) (I) - V^(2) (I)) dl] dx 
As v^(^'(x) > v^(^)(x), we conclude that 



1 

exp 

X Jo 



- r(v/'^^(i)-r<^'(i))df 



dx 



> / exp 

Jo 



V/(i)(x)-v/<2)(^)) 



dx 



dx . 



e Jo 

In particular, pi (x) > x. 

Now, assume that pi(x) — xi > for all x G dS"^K (Note that dS^~^ is an union of a finite number 
of n — 2 dimensional simplexes, where by the principle of induction we assume the result valid.) Finally, 
we use the maximum principle for subharmonic functions to conclude that the minimum cannot be in the 
interior of the simplex (Courant and Hilbert, 1989). Therefore pi (x) >xi for all x G 5"^'. D 

D 

6. The Replicator Dynamics 

In Section 4, we proved that, when genetic drift and selection balance, then there is a special timescale 
such that the evolution of an infinite population can be described by a parabolic partial differential equation. 
Nevertheless, in applications one is usually interested in large but finite populations. In this case, an exact 
limit is not taken, and (18) can be taken as an approximation of this evolution. We shall discuss this further 
in the conclusions, but we observe that this equation might be a good approximation even when balance is 
not exact, i.e., when v and jU are close but not equal to one. This could typically lead to an equation with 
K being either quite large or small. In the former case, a regular expansion in K shows that the evolutions 
is governed by (17). On the other hand, in the latter case, one expects that the much simpler transport 
equation (19) will be a good approximation for the evolution. Indeed, in this section we show that (1') 
can be uniformly approximated by (19) in proper compact subsets of the simplex, and over a time interval 
shorter than K^ ' . 

We start in subsection 6.1 showing that the equation (19) is formally equivalent to the replicator system. 
Afterwards, in subsection 6.2, we answer what we believe to be an important question: what exactly is 
the replicator equation modelling? In particular, we will show, using a simple argument, that the replicator 
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equation does not model the evolution of the expected value (of a given trait) in the population, but the 
evolution of the most common trait (the mode of the probability distribution). Finally, we show, in subsec- 
tion 6.3, that the replicator ordinary differential equation is a good approximation for the initial dynamics 
of the Wright-Fisher process, when K is small. As in Section 5, we shall assume that the fitness functions 
are smooth. 

6.1. The replicator ODE and PDE. We shall now study in more detail the equation (19), which has a 
close connection with the replicator dynamics as shown below: 

Theorem 9. Let 4>r(x) the flow map of 

dx 

(30) d7^^(^^'^^- 

and let 

Q{x,t) = - [\v-a){<i>,^,{x))ds. 
Jo 

then the solution to (19) with a C initial condition po is given by 

(31) p{x,t)=cQ(^''^po{<P-t{x)). 

Proof. Clearly Q{x,0) ~ 0, and <I>o(x) ~ x. Hence the initial condition is satisfied. 
Let 

i?(z,f):=e-e(*'W'')/,(f,4>,(z)) 

On one hand, (31) shows that 



Therefore 



R{z,t)=poiz) 



On the other hand, one can compute 

— (z,f)=e-e(4>rW^')(^^^ + V/?.x+(V.i2)/?)=e-2(*'W.0(5,p + V.(;?i2)). 

We then conclude that p{x, t) is the solution of equation (19). D 

D 

6.2. Peak and average dynamics. We start by showing that the long term dynamics of the average in 
the Wright-Fisher process, even in the thermodynamical limit, is not governed by the replicator equation. 
Consider for example, a population of n types, evolving according to the replicator-diffusion equation with 
fitness functions given by y/*-'' : 5"^' -^ M. 

From the fact that the final state of the the replicator-diffusion equation is given by equation (27), the 
coefficients 7r,[p'], i ~ !,...,« can be calculated in two ways: 

Pi{x)p\x)dx = / pi{x)p°°{x)dx = 7ti[p^] = / Xip°°(x)dx =: {p°°)i . 



Therefore the average of the probability distribution will converge to a certain point of the simplex de- 
pending on the initial condition. This is completely different from the replicator dynamics, as its solution 
converges to a single attractor, periodic orbits, chaotic attractors, etc (Hofbauer and Sigmund, 1998). 

Now, we show that the probability distribution concentrates in the ESS; this shows that the peak will 
behave in manner similar to the solutions of the replicator dynamics. 

Recall that (Hofbauer and Sigmund, 1998), that an ESS that lies in interior of 5"^' must be a global 
attractor of the replicator equation (30). We have then the following result 

Theorem 10. Assume p^{x) is smooth, suppp^(x) C int5"^ and assume that (30) has a unique point x* 
such that for any initial condition x(0) G int5"^', limr^oox(f) = x*. Then the solution of equation (19) is 
such that 

\\mp{x,t) ^8x*. 
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Proof. Assume, initially, that x* G int5"^'. Since x* is a globally stable equilibrium for interior initial 
points, for any given 5 > 0, we can find T > 0, such that, for f > T, we have that for any proper compact 
subset /Tc 5"" ^ 

^,{K)(ZB8{x*). 
Let V^(x) be a continuous function with support contained in K. Then, for t > T,we have that 

p{x,t)\i/{x)dx~ / /?(x,f)v/(x)dx. 
But, let e > be given. Since xj/ is continuous, possibly with a smaller 5 > 0, we must have 

(32) \i/{x*)-e< [ /?(x,f)v/(x)dx< v/(x*) + e, 

Jbs{x*) 

Now take {8k, £k) i such that (32) is satisfied. This yields a sequence of times T^ such that T^-^ oo and 

lim / p(x,ri)v/(x)dx=i/(x*). 

Since ^s{K) C ^t{K), for s> t, the claim follows. 

For the case x* e dS"^\ the result follows from similar arguments, replacing Bs{x*) by Bs{x*) n5""'. 

D 
D 

6.3. Asymptotic approximation. Let 

0<K-<L 

If we perform a regular asymptotic expansion, i.e., if we write p,^ « p'-^^' + Kp^-^^ H , then we find , for 

times f <C ic^ ' , that the leading order dynamics is given by 

(33) drp + V-{pQ.)=0. 

Theorem 11. Assume that the fitness are C (5"^') functions, and that the initial condition p is also 
C {S"^^). Let r^ be the regular part of the solution of (1), with K" > 0. Then po is C (5"^'), and satisfies 
the conservation law (23). Moreover, //" V • i2 > 0, then given K and K positive, there exits a C such that, 
fort <C Ck^ , we have 

\\r^i-,t)-po{;t)\\^<CK 

and 

\\d^po{;t)\\^<K 

Thus pQ is the leading order asymptotic approximation to r^, for f <C K^^C. 

Proof. The statements about po follows straightforward by obtaining the solution by the method of char- 
acteristics. 

Let Wk = r^ — PQ. Then Wk satisfies 

with null initial condition, where 

n—\ n—\ 

g(i{x,t) = ^ df {xipo) - Y, ^ij (xiXjPo) . 

Notice that, because of the assumptions on po, we have that go is uniformly bounded in time. 

The solution for such a problem is given by Duhammel principle. Let S{t,to) be associated solution 
operator. We have that 



K f 

WK{->i,t)^- / S{t,s)g(i{s,x)As. 
2 Jo 



By the maximum principle applied to the semigroup 5(f2,fi), we have that ||5(f,s)^o('Sj-*^)|| ^ M^, and by 
the uniform bound on g^, we have that there exists a constant M such that Ms < M. Thus, we find that 

\\S{t,s)go{s,x)\\^<M. 
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Hence 

\w^{y.,t)\<Kt-. 

Therefore, taking C = 2M"\ we find, for f < CiC"', that: 

||H';c(f,-)ll»<l- 

D 
D 

Remark 8. If the condition on V ■ D, is not satisfied, a similar proof shows that ift<ti — log()c) then the 
same conclusion holds. Notice also that this condition is satisfied if the replicator has a globally stable 
equilibrium in the interior ofS"^. 

Theorem 10 shows that, for sufficient large time, the support of the solution of the replicator PDE, 
equation (19), will be concentrated in sufficient small neighborhoods of x*. In particular, this will be true 
for the maximum. For the replicator-diffusion equation (1) this cannot be valid for any value of K" > (as 
it was proved in theorem 7); however, for strong selection, the initial dynamics given by the replicator- 
diffusion equation is similar to the one given by the replicator ODE. This is justified by the following 
result: 

Theorem 12. Let x^, G S"^ be the unique ESS for the replicator dynamics (30). Let t* be the approximation 
time given by Theorem IL Let £ >0 be given, and assume that there exists a t** <C t* such that the support 
of pQ is contanined in the ball of radius £ around x*. Then, there exists a constant C > such that for 
t** <t < t* we have: 



Be(x.) 



w(x,f)dx-l 



<Ck. 



Proof. By assumption, for t > t**, we have 



Be(x.) 



/?o(x,f)dx= 1 



Let us write Pk — qK + ''K, where r^ is the regular part of /j^-, and sing supp((7^) C dS" ^ By Theorem 11, 
since f <t*, there exists C' > such that ||r^(-,f) -/5o(-,f)ll~ < Ck. But then 



-C'K<rk{x,t)-po{x,t) <C'k , 



and 



Finally 



-Ck< / r^dx-1 <Ck . 

JBeix*) 



/ Pk 

JBeix") 



dx-1 



<Ck. 



a 

D 



7. Numerical results 

We show, in this section, numerical results for two variants of the Rock-Scissor-Paper game (Hofbauer 
and Sigmund, 1998); i.e., fitness are identified with the pay-off from game theory. In subsection 7.1, we 
study the evolution of the discrete evolution numerically in time, and show that the peak of distribution 
behaves accordingly to the replicator equation while the average value of the same distribution converges 
to a point which is not the ESS. In subsection 7.2 we obtain explicitly the fixation probability of a given 
type for the symmetric Rock-Scissor-Paper game. A full animation is available in the website indicated in 
the caption of figure 2. 
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-Sci55ur-Paper game, time = 3 



.: 




.: 










Rock-SciBsor- Paper game, time = 54 






Figure 2. Solution for short times (1,3,6,10,15,21,28,35,44,54,65,77) of the Wright- 
Fisher evolution for a population of 150 individuals of two given types, with fitness given 
by equations (34) and (35) for a distribution initially concentrated in the interior non- 
stationary point j^(70,70, 10). The value of the distribution P{x,y,t) is in logarithmic 
scale. Note that the cyan spot, marking the interior peak of the probability distribution 
rotates and converges to the ESS (5,5,3) (along characteristics of the PDE or, equiva- 
lently, the trajectories of the replicator dynamics). At the same time, the green spot marks 
the mean value of the probability distribution and also rotates initially. After a long time, 
it moves toward its final position, given by x°° := (ci [/?'], C2 [/?'], 1 — ci\p^] —Q [/?']) ~ 
(0.331,0.227,0.442). For a full animation, also for different population sizes A^, see 

http : //dl . dropbox . com/u/11325424/WFsim/RSPFinal .html 



7.1. Forward equation. We use evolutionary game theory (Smith, 1982; Hofbauer and Sigmund, 1998) 
to define the fitness function. More precisely, we define a pay-off matrix M — (Mij). .^j ^^ such that Mij 
is the gain (in fitness) of the / type against the j type. The fitness of the / type in a population at state x is 



(34) 



V/«(x) = £m,-,-^,-(Mx)^. 
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a. 8 



0.4 



Figure 3. Fixation probability of the third type, in a Rock-Scissor-Paper game. This is 
the numerical solution of the stationary state of the equation (28), simulated by a Wright- 
Fisher process with A^ = 150 and pay-off matrix ( [[20, 0, 40] , [40, 20, 0] , [0, 40, 20]] ) . Note 
that higher values of the fixation probability "rotates" around the center of the simplex 
(the stationary state of the replicator dynamics). 



The replicator dynamics is given by the system of differential equation i,- — x,(v/'''(x) — V/(x)), where 
\j/{x) = X ■ Mx. 

We consider in Figure 2 the evolution of a discrete population of N = 150 individuals with the pay-off 
matrix given by 



(35) 



M: 



30 81 29 

6 30 104 

106 4 30 



This is know as the generalized Rock-Scissor-Paper game and presents an evolutionary stable state (ESS) 
{x* ,y* ,z*) ~ (j; jj 5)- Furthermore, the flow of the replicator dynamics converges in spirals to the ESS. 
The vertices as well as (x*,3'*,z*) are equilibrium points for the continuum dynamics. See Hofbauer and 
Sigmund (1998) for the choice of values of the matrix M. 

Note that the peak moves in inward spirals around the central equilibrium, following the trajectories of 
the replicator dynamics, while all the mass diffuses to the boundary. 

The green spot indicates the average value for x and y\ at first it moves in spirals close to the trajectories 
of the replicator dynamics. After a time depending on the value of A^ it starts to move in the direction of its 



final point (x" 



°) = (tTi [/?'], 712 [p^], 7^3 [p']). This point can be calculated using equation (27) and the 



n — T) independent conservation laws. Effectively, let x, denote a given vertex of the simplex 

lim(x,)(r) — lim / Xip{t^x)Ax — ni\p ] , 

where %i \p^\ is the fixation probability of type i associated to the initial condition p^. 

1.2. Backward equation and the decay of the interior L^-norm. The stationary state of the backwaid 
equation (28) represents the fixation of probability of a given type. This type is specified by the associated 
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10 20 30 

time-steps 



40 



Figure 4. Evolution of the probability mass, for the Rock-Scissor-Paper game given 
by matrix (36) and with initial condition concentrated in the ESS, p' = 8,[ ly The red 

line indicates the mas {L^ -norm) in the interior of the simplex; the blue line, the mass in 
the interior any of the faces, and the black line, the mass in any of the vertices. 



boundary conditions. Let us consider, as an example, that n ■ 
Paper game defined by the matrix 



: 3, the evolution is given by the Rock-Scissor- 



(36) 



40 20> 
M = I 20 40 
i40 20 



and we study the fixation probability of the third type. An exact solution is difficult to obtain, as it would 
be necessary to solve an hierarchy of equations, each solution representing boundary conditions of a larger 
set; however, a numerical solution is extremely easy to compute, as the Wright-Fisher process is a natural 
discretization of the (forward as well as the) backward equation (cf. theorem 5). This is probably compu- 
tationally inefficient, and different processes can be compatible with the same limit equations. See figure 3 
for an illustration. 

In figure 4, we plot the L' norm in the interior of the simplex and all subsimplexes, showing that that 
the probabiUty mass flows from the simplex 5"^' to the faces (which are equivalent to the simplexes S"^^); 
the solution behaves on the faces as the solution of the replicator-diffusion problem with one dimension 
less. The probability flows to the "faces of the faces", i.e., to simplexes S"^^ until it reaches the absorbing 
state e, (simplexes S^) for / = 1 , . . . , n. We may think in a stochastic process reaching and sticking to the 
faces of the simplex until they reach their final spot, the vertices; note that accoding to figure 4 there is a 
non-negligible probability of simultaneous extinctions. 
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8. Conclusions 



We present a derivation of continuous limits of discrete Markov chain evolutionary models, that are 
frequency-dependent extensions the classical Wright-Fisher model, through pure analytical techniques. 
The derivation presented pays close attention to the variety of possible time scalings possible as related to 
the selection and population size that are measured by two parameters /i, v > 1. The balance of diffusion 
and selection (/i = v = 1 in our terminology) can be seen as slight extension of the results in (Ethier and 
Kurtz, 1986, Chapter 10) using analytical methods instead of probabilistic arguments, and that favours the 
forward Fokker-Planck equation instead of the backward. In this sense, from a mathematical point of view, 
the weak formulation presented for the forward equation seems to be new, in particular for the minimal 
assumptions on the fitness functions. The case 2> jj. > v — I yields a hyperbolic equation that is the PDE 
version of the replicator equation. An apparently similar result can be found in (Ethier and Kurtz, 1986, 
Chapter 11), which would correspond formally to take Af = 1 and N ^ I, without using explicit scaling 
between these two variables. 

With some additional regularity assumptions on the fitnesses functions, we can show that (18) is equiva- 
lent to (1') together with the conservation laws (23). In particular, this allows to caracterise the behaviour of 
p on the lower dimensional subsimplices of S"^ ^ . This can be used to obtain equations for the probability 
of extiction among other information. 

The results here are also related to results in Champagnat et al. (2006, 2008), where the idea that the 
underling scaling influences the macroscopic model was already present, although in a less explicit way 
than here. Nevertheless, the accelareted birth-death regime in Champagnat et al. (2008) can be seen as 
a counterpart to our scaling of At and A^. On the other hand, the scaling for the fitness are taken as fixed 
(corresponding to our v = 1 in our terminology), and this explains why they do not obtain the pure diffusive 
limit in the large population regime. Notice also that the large population regime taken there seems to 
anihilate any stochastic effects coming from births and deaths, and that the stochastic effects in this limit 
are due only to the mutation process. 

However, as pointed out above, as we allow more flexibility in the scaling laws, we are able to highlight 
any of these two factors independently; more precisely, for certain choices of the scaling in the fitnesses 
functions (namely, the exponent v), their influence in the dynamics goes to zero so fast that the limit 
model is purely diffusive. On the other hand, if we grow the population size fast enough (i.e., jU > 1) then 
we higlight the determinist evolution, providing a direct way to compare the replicator equation with the 
Wright-Fisher process (or, for that matter, also with the Moran process, but, naturally, in a different time 
scale). To the best of our knowledge, this explicit comparision is new. See also Fournier and Meleard 
(2004) for a similar approach. 

The use of ordinary differential equations in population dynamics is widespread. However, as they are 
valid only for infinite populations, and real populations are always finite, the precise justification of its 
use and the precise meaning of its solution is seldom made clear. In this paper, we showed, in a limited 
framework, but expanding results from previous works (Chalub and Souza, 2009a,b), that ODEs can be 
justifiably used to model the evolution of a population. However, the validity of the modeling is necessarily 
limited in time (increasing with the population size), and the solution of the differential equation models 
the most probable state of the system (therefore, the differential equation would give answers compatible 
with the maximum likelihood method, but not necessarily compatible with other estimators). 

One of the central issues of the present work is to discuss the possibility of using diffusive approximation 
for large, but finite, A^. On the other hand, a major challenge to any one interested to use the replicator- 
diffusion equation to fit experimental data is the value of the K. However, we could give some hints in 
its meaning. Assume t is finite and A^ is large enough. In this case, both the replicator and the replicator- 
diffusion equations are essentially indistinguishable. Therefore, K should be small. In a sense, K^^ is a 
possible definition of effective population size (see also Etheridge (201 1) for alternative definitions). 

On the other hand, what in this work is called scaling consists in the triple {ij.,v,k) and these are 
essentially independent parameters. Assume v = 1 and /x is slightly superior to 1. In this case, the thermo- 
dynamical limit is the replicator equation; however, for finite / and A^ large enough, one might expect that 
the solution cannot be far from the solution of the replicator-diffusion equation; we conclude again that K 
should be small. This argument also apply to the case /i = 1 and v — 1 ^ 1. 
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Finally, we would say that in the case where the population is large, time is limited, and natural selection 
and random genetic drift are essentially balanced, it is natural to expect a small value of K. If natural 
selection is a dominant effect, then we also expect a small value of K. Only when the population is small 
or times are long in the evolutionary scale, we would expect order 1 values of K. 

We are currently applying a similar technique to epidemiological models; in this case it is necessary to 
impose boundary conditions in part of the boundary (as an homogeneous population of infected individual 
is not stationary, as infected individuals become, with time, removed or even susceptible) and it is impossi- 
ble to impose boundary conditions in part of the boundary (a population of susceptible remains in this state 
for ever). Early results were already published in Chalub and Souza (2011). The same problem, regarding 
the imposition of boundary conditions is true if we include mutations in the Moran or Wright-Fisher model. 
This is work in progress. 
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Appendix A. Integral (10) 

Let F be an (n — 1 ) X (n — 1 ) matrix and x e M'j_ a vector such that Fa =x^ +x^^ and Fij = x^\ Then 
detF = Y/X, where Y = L"=i x,-, X = n"=i ^^ 

Now, let G be a positively defined matrix, and ^ its associated quadratic form. Then 



T,-T,-e-2'^(^'^)dT ^ da, J e-z-^^'^)dT = (2;r)^5G,, 



2J'"' "^"j ' ' ^'^VMg 



, "-1 



2(detG)J 



(-l)'+^detG' 



'J) 



where G'^'-'' is the matrix obtained from G when we remove the /-row and the y-column. 
Imposing G = F, / = j and F = 1 we find 



/ 



T2e-z*(^-^)dT= (2;r)^Xi(-l)2'§^M^ 



. n—l 3 X; { i — X; ] , ,1—1 / — 

{2n)^X^ "■ " = {2n)^VXxi{l -x; 



Now, let us assume / =^ j, and consider the matrix F'-'-''. The {j,j) entry of the original F matrix was 
removed; therefore, the y'-row is transformed in a row made only of x^' entries; equally the (/, /) entry was 
removed, therefore there is a column made only of x^^. Subtract x^' to every other row; the said column 
is transformed such that it has now only zeros, except for one entry, equal to x^ ' . Finally, 

X ' 



detF("' = ^^ n xk] =(-1) 

^" \ke{l,...,n-l}\{i,j} 



and this finishes the proof. 



Appendix B. Weak formulation in time 



In order to obtain a truly weak formulation, without any requirement upon the regularity of p, we 
observe that the equation above is valid for any time t/^ — to + kAt. Hence, if we also let T = {m+l )At in 



32 



FABIO A. C. C. CHALUB AND MAX O. SOUZA 



the equation above, and sum over k, we obtain that 



L E iPNi'>^,tk+l)~PNi-!l;tk))8ix,tk) 



1 

2N 



n-l 



E E PN{^,tk){ E ^'(5o-^>)^;y-^(x,f^) 



1^=0 xes"-' 



Y^iAty Y. p(^''^) 



'n-l 



k=^ XGS- 



>1 



dx. 



On summing by parts the left hand side, we obtain 



m— 1 



Y. Y. PN{->i,tk){g{x,tk+i)-g{x,tk)) 



^=OxGS"-' 



Y PNix,tQ)g{x,to)dx+ Y, /'A'(x,r)^(x,r) 



xGS" 



xes'L 



1 m I n—\ \ 

^E E /'A'(x,f-t) E -*^'(^'7"-*^/)^o'^(^'^*) 



-^(Af)^ Y Pi^jk) 

k=0 xGS"-l 



1-1 



YxiU^^H^)-¥{^))di8{xA) 

,.7=1 
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